The main transition in the Pink membrane model: finite-size scaling and 

the influence of surface roughness 
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We consider the main transition in single-component membranes using computer simulations of 
the Pink model [D. Pink et al, Biochemistry 19, 349 (1980)]. We first show that the accepted 
parameters of the Pink model yield a main transition temperature that is systematically below 
experimental values. This resolves an issue that was first pointed out by Corvera and co-workers 
[Phys. Rev. E 47, 696 (1993)]. In order to yield the correct transition temperature, the strength 
of the van der Waals coupling in the Pink model must be increased; by using finite-size scaling, a 
set of optimal values is proposed. We also provide finite-size scaling evidence that the Pink model 
belongs to the universality class of the two-dimensional Ising model. This finding holds irrespective 
of the number of conformational states. Finally, we address the main transition in the presence 
of quenched disorder, which may arise in situations where the membrane is deposited on a rough 
support. In this case, we observe a stable multi-domain structure of gel and fluid domains, and the 
absence of a sharp transition in the thermodynamic limit. 

PACS numbers: 87.16. D-, 87.14.Cc, 64.70.-p, 82.20.Wt 



I. INTRODUCTION 

Lipid membrane bilayers are abundant in nature and 
to understand their properties is of paramount impor- 
tance [IH1] • One aspect that has received much attention 
are collective phenomena (phase transitions) taking place 
in these systems. Among the different phase transitions 
that can occur [143 > the main phase transition ispre- 
sumably the most important and well studied one 0, [t| . 
This transition, typically driven by the temperature T, 
is between a "gel" and a "fluid" phase. At low T, the bi- 
layer is in the gel phase (characterized by nematic chain 
order of the lipid tails), while at high T the bilayer as- 
sumes the fluid phase (characterized by the absence of 
nematic chain order). 

Computer simulations have become a well established 
tool to model the main transition. The challenge in sim- 
ulations is to strike a balance between the level of detail 
to include, and the time and length scale one wishes to 
address [1 01 ] . Since collective phenomena involve many 
molecules and entail large length scales it is clear that, 
in order to describe the main transition, a significantly 
coarse-grained particle model is crucial. Strictly speak- 
ing, one needs to address the thermodynamic limit (in- 
finite particle number) since only there phase transition 
properties become properly defined. Indeed, the need for 
coarse grained modeling of lipid bilayers is well recog- 
nized [1M3. 

An early and highly successful coarse grained approach 
to study the main transition has been the particle model 
introduced by David Pink and co-workers [Lll - tla ] . In this 
model, the so-called Pink model, only the orientational 
degrees of freedom of the hydrophobic lipid tails are in- 
cluded, while the positional degrees of freedom of the hy- 
drophilic heads are disregarded. This model, due to its 
simplicity, allows for the investigation of very large sys- 
tems, and the nature of the main transition can be probed 



in great detail. Indeed, key features of the main transi- 
tion in the Pink model compare well to experiments [l7| ■ 

However, despite the great success the Pink model has 
enjoyed, there remain some open questions. In Ref. [l8l . 
it was noted that the Pink model at the experimentally 
determined transition temperature does not undergo any 
transition. While in systems of finite size there were indi- 
cations of a transition, these vanished in larger systems. 
This raises the question as to why no transition could 
be detected. The aim of this paper is to resolve this is- 
sue. As it turns out, to properly model the main transi- 
tion, a finite-size scaling study is essential. Computer 
simulations inevitably deal with only a finite number 
of particles, and their output will depend on the num- 
ber of particles used, especially near phase transitions. 
Finite-size scaling provides the framework to systemati- 
cally extrapolate simulation data to the thermodynamic 
limit. To date, finite-size scaling studies of the Pink 
model are scarce, with Ref. [H being a notable excep- 
tion. The present paper aims to fill this gap. Our main 
finding is that, in order to observe the main transition 
in the Pink model at experimentally relevant tempera- 
tures, one of the model parameters needs to be adjusted. 
This follows quite naturally when one realizes that the 
universality class of the Pink model is just the one of 
the two-dimensional (2D) Ising model [la]. As we will 
show for three lipid species, the "standard" Pink model 
parameters yield a critical temperature distinctly below 
the experimental main transition temperature. Conse- 
quently, a "re-tuning" of the standard Pink parameters 
is urgently needed. 

As an application, we also address the fate of the main 
transition in the presence of quenched (immobilized) im- 
purities using the Pink model. The experimental motiva- 
tion to do so is that this situation may resemble that of 
a membrane supported on a rough substrate. In binary 
lipid mixtures, the effect of such impurities on lateral 
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phase separation has recently attracted much attention 
[19l - l24| . In this paper, we present simulation results for 
the corresponding scenario in a single component bilayer 
undergoing the main transition. Within the framework 
of the Pink model, we find that quenched impurities pre- 
vent the main transition from taking place, already at 
low impurity concentrations. Instead of the formation of 
macroscopic gel and fluid domains, we now obtain a sta- 
ble multi-domain structure, which strikingly resembles 
experimental results. The theoretical justification is that 
the impurities induce a change in universality toward the 
2D random-field Ising class. As is well known, the latter 
does not support an order-disorder phase transition in 
the thermodynamic limit [25|-[28j|. 



II. THE PINK MODEL 

In the Pink model, the lipid bilayer is assumed to con- 
sist of two independent monolayers. Each monolayer is 
represented by a triangular 2D lattice consisting of N 
sites, and each lattice site contains a single lipid chain. 
Each lipid molecule is comprised of two independent hy- 
drophobic acyl chains and a hydrophilic polar head. The 
polar heads are translationally frozen to the lattice, and 
no particular structure for the polar head groups is as- 
sumed. The only degrees of freedom included in the Pink 
model are the acyl chain conformations. These are not 
simulated directly (i.e. one does not explicitly model the 
carbon atoms) but are captured in a coarse-grained fash- 
ion whereby the chain conformations arc grouped into 
a = l,...,q discrete states. The original Pink model 
uses q = 10, but we will consider different values also. 
These states include the ground state (a = 1), eight low- 
energy excitations (a = 2, . . . , q — 1), while all remaining 
conformations are grouped into a single disordered state 
(a = q). Each state a is characterized by three coarse- 
graining parameters, namely an internal energy E a , a 
cross-sectional area A a , and a degeneracy D a counting 
the number of chain conformations with energy E a and 
area A a . 



A. coarse graining parameters 

To determine the coarse graining parameters, we as- 
sume that a single acyl chain consists of i = 1, . . . , M 
carbon atoms, thereby containing M — 1 carbon-carbon 
bonds, and that bonds are either in a trans or gauche 
configuration. The trans configuration yields the low- 
est energy, while the gauche configuration has a slightly 
higher energy. The energy difference between the trans 
and gauche configuration is denoted T (Table HJ. To un- 
derstand the difference in geometry between trans and 
gauche bonds consider a chain segment of four consecu- 
tive carbon atoms. The positions of the first three atoms 
define a two-dimensional plane. In the trans configura- 
tion, the fourth atom remains in the plane, while in the 




5) / 




©X® 








f) ( 


7 






(b) 



FIG. 1. Typical chain conformations of the Pink model with 
M = 7, showing (numbered) carbon atoms placed on the 
nodes of a hexagonal lattice. The atom connected to the 
head group is labeled i — 1 but for clarity the head group is 
only drawn for the ground state. The z-direction indicates 
the bilayer normal, while the vertical double-arrows indicate 
the projected length, (a) The ground state a = 1, consisting 
of only trans bonds, (b) The two conformations that consti- 
tute the first excited state a — 2 containing one gauche bond 
(marked with a cross), (c) Conformation belonging to the 
second excited state a — 3. The internal energy is the same 
as in (b) but the projected length is shorter (the other a = 3 
conformation has the gauche bond between atoms 3 — 4) . 



gauche configuration, it leaves the plane, and it can do so 
inward or outward. Thus, each gauche bond is two-fold 
degenerate. In the Pink model, it is assumed that each 
2n-th gauche bond takes the chain back to the original 
plane, and so the gauche degeneracy is given by 



q 2 ce il( n /2) 



(1) 



where n denotes the total number of gauche bonds in the 
chain, and where the function ceil means "rounding-up" 
to the nearest integer. 

It is convenient to mathematically represent the chain 
conformations on a hexagonal lattice with next-nearest 
neighbor distance 2a. We emphasize that this lattice 
is merely an aid to identify the low energy chain con- 
formations which are needed to set the coarse-graining 
parameters: it should not be confused with the trian- 
gular simulation lattice on which the Pink Hamiltonian 
will eventually be defined. The carbon atoms are placed 
on the nodes of the hexagonal lattice following certain 
rules, and nearest-neighbor connections between atoms 
represent carbon-carbon bonds. The ground state a = 1 
corresponds to the chain conformation that is maximally 
stretched [Fig.QTa)]. Note that, in the ground state, the 
atoms are alternatingly placed on the left and right lattice 
node, yielding a characteristic "zig-zag" pattern. The 
ground state by definition contains only trans bonds, its 
internal energy is set to zero as a reference E\ = 0, and it 
is obviously non-degenerate D\ = 1. The cross-sectional 
area of the ground state has experimentally been deter- 
mined as Ai = 20 .4 A 2 [TJ]. We also introduce the pro- 
jected length I of the conformation, defined as the dif- 
ference in z-coordinate between the carbon atom closest 
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to the head group (i = 1) and the one furthest away 
(i = M), with the z-direction as indicated in the figure. 
For the ground state, it follows that l\ = (M — l)a. 

The eight low-energy excitations (a = 2, ... ,9) are 
obtained by systematically incorporating gauche bonds. 
The effect of such a bond is to disrupt the "zig-zag" pat- 
tern of the ground state. That is, one no longer places the 
atoms altcrnatingly on left and right nodes, but also al- 
lows for "excursions" whereby for two consecutive atoms 
the same direction is chosen. Each such excursion cor- 
responds to a gauche bond, and has energy cost T. The 
gauche bonds are introduced according to the following 
rules: (1) The two bonds in the chain closest to the 
head group must always be in the trans configuration. 
In Fig. [U these correspond to the bonds between atoms 
1 — 2 and 2 — 3. (2) At most three gauche bonds are 
allowed, and each time such a bond is included there is 
an energy cost T. (3) The projected chain length I must 
obey h — I < 3a. (4) The acyl chain cannot fold back 
onto itself. In the coordinate system of Fig. [TJ this means 
that the z-coordinates of the atoms must obey 2^+1 > z,. 

Following these rules, we show in Fig. [TJb) the chain 
conformations (i) and (ii) that form the first excited state 
a = 2. In (i), a single gauche bond is placed at the 
very chain end, while in (ii) it is placed at the second- 
last position. One immediately sees that both confor- 
mations have the same energy E 2 = T, and the same 
projected length l 2 = {AI — 2)a. To compute the cross- 
sectional area, one assumes volume conservation for the 
lipid chains: A a l a = Ail±. Hence, from the (known) 
ground state values, the cross-sectional Ml CM. of the ex- 
cited state follows. Note that, by placing the gauche 
bond at the third- last position [Fig. (He)], a shorter 
projected length is obtained, and so conformation (c) 
does not belong to the first excited state (even though 
it has the same energy). The total degeneracy of the 
first exited state D2 = 4, which is the total number 
of conformations, multiplied by the gauche degeneracy 
of Eq. ([T]). The coarse- graining parameters of the re- 
maining excited states can be found analogously, and 
arc listed for completeness in Table |TJ Finally, for the 
completely disordered state a = q = 10, one assumes 
E w = (0.42M - 3.94) x 1CT 13 erg, A 10 = 34 A 2 , and 
degeneracy D w = 6 x 3 M ~ 6 , which have their origins in 
experimental considerations (l6j | . 



B. Pink model Hamiltonian 

Having specified the coarse-graining parameters, the 
Hamiltonian of the Pink model can be written as [3fJ 

'Hpink = Hq + Hvdw + Hp ■ (2) 

The first term is the total internal energy of the acyl 
chains %q = 53 i=1 E s u\, with the sum over all N sites 
of the triangular lattice, and s{i) £ {1, . . . , q} the confor- 
mational state at the i-th lattice site. The second term 
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TABLE I. The coarse graining parameters used to describe 
the acyl chain conformations in the q — 10 Pink model p^-[Tf3. 
[2S|. For each state conformation a, we list the internal energy 
E a , the projected length l a , and the degeneracy D a . The 
energy of a single gauche bond equals F = 0.45 x 10~ 13 erg, 
while M denotes the number of carbon atoms in the chain. 

represents the anisotropic van der Waals interaction be- 
tween adjacent acyl chains "Hvdw = -^oE(! j\ ^s(i) ^s(j) 1 
with Jo the van der Waals coupling constant, and 
a sum over all 3./V nearest-neighboring sites on the tri- 
angular lattice. The precise value of Jo depends on the 
chain l ength and explicit expressions are provided else- 
where [l4l [l5l.[3l|. However, it has been noted that these 
parameters do not always yield a main transition at the 
expected temperature [18|, and so we will also propose 
our own values later on. The (dimensionless) variables 
I a measure nematic chain order, and can be expressed in 
terms of the cross-sectional areas 0, [l8[ 

(HH)(*r ■ < 3 > 

where ujiq = 0.4 for the disordered state a = 10, and 
uj a = 1 otherwise. 

The last term in the Hamiltonian accounts for the in- 
teraction between the hydrophilic polar head groups and 
between them and water and also steric interactions from 
both head groups and the lipid chains. Although it is 
possible to consider a more realistic pairwise interaction 
between the headgroups [30| , this interaction can be ap- 
proximated with a simple pressure term Hp = HA, where 
n is an effective lateral pressure acting on the lipid chains 
in the bilayer membrane, and A the total cross-sectional 
area occupied by the lipids chains 

N 

a = Y,M) ■ ( 4 ) 

i=l 

III. MONTE CARLO METHODS 

To study the phase behavior of the Pink model, we use 
the Monte Carlo (MC) simulation method. We mostly 
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use triangular lattices of size N = L x L with periodic 
boundary conditions. The principal MC move consists 
of randomly picking one of the lattice sites, read-out the 
conformational state a of that site, and propose a new 
state /3 drawn randomly from the set of q possible states. 
The new configuration is accepted with the Metropolis 
criterion 



(a) 



Dp (An 

1, — exp ~~r~Pp 

D a V k B T 



(5) 



where D denotes the state degeneracy, AH the energy 
difference between initial and final configuration as given 
by Eq. k B the Boltzmann constant, and T the tem- 
perature. The degeneracy compensates for the fact that 
some of the states have a much larger entropy, and should 
therefore appear more often in the ensemble average. 

By virtue of the MC move, the total projected area 
A given by Eq. (0} fluctuates during the course of the 
simulation. In fact, A plays the role of order parame- 
ter since it changes abruptly at the main phase transi- 
tion. Hence, it is instructive to measure the distribution 
P(A\T, II), defined as the probability to observe a con- 
figuration with projected area A. The distribution de- 
pends on the imposed temperature and pressure, as well 
as on the linear extension L of the triangular simula- 
tion lattice. At the main transition, P(A\T,U) assumes 
a characteristic bimodal shape, from which a number of 
important phase properties are obtained (explicit exam- 
ples are provided in the next section) . We note that even 
with very long simulation runs, distributions P(A\T,TT) 
of high statistical quality are difficult to obtain, espe- 
cially in the vicinity of the main transition. The reason 
is related to free energy barriers that arise from the for- 
mation of interfaces [32h34| . To overcome this problem, 
we combine our MC simulations with a biased sampling 
scheme called successive umbrella sampling [35| ; the lat- 
ter ensures that P(A\T,H) is sampled accurately over 
the entire (specified) range in A of interest. A final in- 
gredient to economize simulation time is the use of his- 
togram reweighting [36j. A single simulation run yields 
P(A\T, IT) at a given temperature T and effective pres- 
sure II; histogram reweighting enables us to extrapolate 
the measured distribution to different values T', II'. For 
example, extrapolations in the pressure are performed 
using P(A\T,W) oc P(A\T, II) e -(n' ~n)A/k B T _ Extrapo _ 
lations in the temperature can be performed analogously, 
but also require storage of the energy histograms; for im- 
plementation details see Ref. H3- 



IV. RESULTS 

A. the "standard" Pink model revisited 

We first consider the main transition in a membrane 
consisting of DPPC lipids to settle a controversy when 
this system is being simulated using the Pink model. The 
acyl chains in DPPC consist of M = 16 carbon atoms, 
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FIG. 2. Simulation results for DPPC obtained using the Pink 
model with "standard" parameters, (a) Probability distribu- 
tion P(A) of the cross-sectional area per molecule. Note that 
we have adopted the convention to plot the average area per 
lipid, A = 2A/N, with A given by Eq. (g]|. At high temper- 
ature, irrespective of the value of II, P(A) is single-peaked 
corresponding to one phase (solid line). At low temperature, 
P(A) becomes double-peaked provided II = IIcoex, indica- 
tive of two-phase coexistence (dotted line), (b) Finite-size 
scaling analysis to locate the critical temperature T c . Plotted 
is the Binder cumulant Ui as a function of temperature T for 
different system sizes L. The intersection of the curves for 
different L yields T c . 



and the experimentally obtained main transition tem- 
perature Tdppc = 314.0 K [3l|. However, simulations 
based on the Pink model could not detect a transition 
at this temperature (l8l | . The latter simulations used the 
"standard" Pink parameters as listed in Table HI van der 
Waals coupling constant J = 0.710 x 10~ 13 erg, and 
pressure n = 30 dyn/cm. Hence the question arises 
as to why no transition could be detected. To answer 
this question we perform additional DPPC simulations 
usin g th e Pink model, with the same parameters as in 
Ref. [13, but over a wider range in temperature and pres- 
sure. The picture that emerges is the following: At high 
temperature the distribution P(A\T, II) is always single- 
peaked (corresponding to one phase) for all value of the 
lateral pressure n. At low temperature, P(A\T,H) is 
doublcd-peaked for a special value of the lateral pres- 
sure, n = ncoEX, corresponding to two-phase coexis- 
tence [Fig. Ufa)]. Here, the left peak reflects the gel- 
phase, the right peak the fluid-phase. The numerical 
criterion to locate Hcqex is to vary n until the fluctua- 
tion (A 2 ) — (A) 2 reaches a maximum [38| . with the ther- 
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TABLE II. Critical point parameters for three lipid species, 
with M the number of carbon atoms in a single chain. We list 
the critical temperature T c and coexisting pressure IIcoex 
obtained in simulations of the Pink model using the "stan- 
dard" value of the van der Waals coupling constant Jo. The 
resulting estimates of T c are to be compared to the experi- 
mental melting temperatures T m : T c clearly underestimates 
T m in all cases. Instead, by using the Pink model with the re- 
tuned values Jo proposed in this work, T c coincides with T m , 
with corresponding critical pressure IIq OEX (coupling con- 
stants in units of 10 -13 erg, temperatures in K, and pressures 
in dyn/cm). 
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FIG. 3. Susceptibility scaling function X-t -7 '" versus tL 1 ^" 
for DPPC obtained using the "standard" Pink model. The 
data for different system sizes strikingly collapse using 2D 
Ising values for the critical exponents. 



mal averages computed as (A m ) = / A m P(A\T, II) dA. 

At the temperature T = T c where the transition from a 
single to doublcd-pcakcd distribution occurs, the system 
becomes critical. To locate the critical temperature a 
finite-size scaling analysis is performed, whereby we plot 
the Binder cumulant U\ = (A 2 )/(|A|) 2 , A = A - (A), 
versus temperature T for different system sizes L. In the 
thermodynamic limit 

r i t < t c , 

lim U 1 = {Ui T = T C1 (6) 
[tt/2 T>T c , 

while in systems of finite size, curves for different L inter- 
sect at T = T c Hi|. In Fig. we show the result 
for DPPC obtained using the "standard" Pink model pa- 
rameters: the data scale as expected, and from the in- 
tersection the critical temperature T c can be accurately 
"read-off" . 

The corresponding estimates of T c as well as the coex- 
istence pressures IIcoex for three lipid species are col- 
lected in Table HO For all lipid species considered, the 
computed critical temperature T c is distinctly below the 
experimental melting temperature T m . In other words: if 
one simulates the Pink model at the experimental melt- 
ing temperature T m , one is always inside the one-phase 
region, where P(A\T, II) is single-peaked! This, appar- 
ently, is the reason why no phase transition could be 
seen in previous studies [l8|]. One possibility to get the 
proper value for the transition temperature, i.e. such that 
T c coincides with T m , is to re-tune the value of Jo- This 
has been done for the three lipid species by systemati- 
cally changing the coupling constant Jo using histogram 
reweighting and finite-size scaling. Our proposed values 
Jg and corresponding pressures n coEX for the three lipid 
species are summarized in Table HT1 

For completeness, we still confirm the universal- 
ity class of the critical point, which for the Pink 
model is expected to be the one of the 2D Ising 



model [Ijl. To this end, we consider the susceptibil- 
ity X = ((A 2 ) - (|A|) 2 ) /(k B TL 2 ) [IH, which diverges 
at the critical point \ cx |£| -7 , t = T/T c — 1, with critical 
exponent 7. In systems of finite size, the divergence is 
rounded, but 7 can still be obtained using the standard 
finite-size scaling procedure of plotting x L~ 1 / u versus 
tL 1 '" [13], where v is the correlation length critical ex- 
ponent. Provided suitable values 7, v, T c are used, data 
for different L collapse. The result for DPPC is shown 
in Fig. 02 where the "standard" parameters of the Pink 
model were used. Indeed, by using the 2D Ising val- 
ues {7 = 7/4, v = 1}, and T c = 291.7 K of Tabic M 
an excellent data collapse is observed (similar good col- 
lapses arc obtained for DMPC and DSPC also). The or- 
der parameter critical exponent has also been measured, 
and the 2D Ising value f3 = 1/8 was confirmed (scaling 
plot not shown). Therefore, even though the Pink model 
is a 10-state model, its critical behavior remains in the 
universality class of the 2D Ising model. This further 
motivates the idea of reducing the q = 10 states in the 
Pink model to an effectively two-state description as is 
frequently done 0, d ■ 



B. modified Pink model with fewer states 

We now consider the effect of lowering the number 
of states in the Pink model. For this purpose, an ap- 
propriate number of intermediate states was removed, 
based on the maximum number of gauche bonds. In the 
"standard" 10-state Pink model at most three gauche 
bonds are allowed. Wc now consider the case where 
at most two gauche bonds are permitted, by removing 
states a = 8, 9 from Table [I] yielding an 8-state model 
(to keep the total number of states constant the degen- 
eracy of the removed states was added to the disordered 
state, but we emphasize that this correction is small). 
We apply our previous finite-size scaling analysis to the 



resulting 8-state model for DPPC, using the "standard" 
value J = 0.710 x 10~ 13 erg. As expected, the crit- 
ical point remains in the universality class of the 2D 
Ising model, but it is "shifted" to T c = 309.4 K and 
IIcoex = 26.0 dyn/cm. Similarly, by allowing at most 
one gauche bond, a 5-state model is obtained. In this 
case, the DPPC critical point is located at T c = 351.5 K 
and IIcoex = 87.7 dyn/cm. 

Therefore, lowering the number of states in the Pink 
model while leaving the other parameters untouched, one 
finds that both the critical temperature and pressure in- 
crease. This trend is consistent with the Pink model 
simulations of Ref. HH for DPPC performed at the exper- 
imental melting temperature T m but with a lower number 
of states. In these simulations, hysteresis loops indicat- 
ing a first-order transition are clearly visible around T m 
for q < 6. Indeed, as our scaling analysis shows, by low- 
ering the number of states q, T c eventually exceeds T m , 
resulting in a genuine phase transition at T m . 

To conclude: lowering the number of states q does not 
affect the universality class of the Pink model, which re- 
mains 2D Ising (provided q > 2, of course). Hence, the 
topology of the phase diagram remains the same, merely 
the critical point gets shifted. Depending on the parame- 
ters used, the main transition in the Pink model is cither 
first-order (T < T c ), or it is 2D Ising critical (T = T c ). 
We do not claim that the main transition as observed 
in experiments necessarily conforms to this scenario (we 
return to this point in Section V). 
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FIG. 4. DSPC simulation results for T — 291 K in the absence 
of quenched disorder, (a) The natural logarithm of P(A) at 
IIcoex and system size L = 50. The barrier AF is related to 
the line tension via Eq. (0. (b) Typical snapshot of the bi- 
layer with the lipids color-coded according to their conforma- 
tional state for a 200 x 300 lattice. The snapshot was taken at 
cross-sectional area A = 54.4 A 2 chosen "between the peaks" 
of P(A). A pronounced coexistence between a single gel and 
fluid domain is observed. 



C. Pink model with quenched disorder 

As a final illustration, we consider the main transition 
in a solid-supported membrane, which has received con- 
siderable attention in experiments [H, H3-IEH . A striking 
feature observed in one of these studies is the formation of 
coexisting gel and fluid domains that do not coalesce with 
time, but instead form a multi-domain structure that is 
stable over hours To understand the stability of 

this structure is not trivial, due to the large amount of 
line interface it contains. Here we attempt to reproduce 
such a multi-domain structure within the framework of 
the Pink model. Our hypothesis is that the solid sup- 
port onto which the membrane is deposited has a certain 
roughness. Since surface roughness is random and time 
independent, it constitutes a form of quenched disorder. 
We assume that this gives rises to regions on the sur- 
face where certain lipid tail conformations are preferred 
over others. We capture this effect in the Pink model by 
randomly labeling a fraction p of the lattice sites as "pin- 
ning sites" . At the pinning sites, the corresponding lipid 
chain is fixed into the ground state conformation. This 
extension is trivially incorporated into our MC simula- 
tions: we simply do not apply the MC move to pinning 
sites. We specialize to DSPC, using the "standard" Pink 
parameters of Tables [T] and [TTJ 

In Fig. Ufa), we show lnP(A\T,IL) at T = 291 K and 



ncoEX = —18.9 dyn/cm in the absence of quenched dis- 
order (p = 0). At this temperature, which is well below 
T c , the main transition is strongly first-order. Conse- 
quently, there is a significant line tension a between gel 
and fluid domains; the latter is related to the free energy 
barrier [13, HI 

AF = k B Tln{P max /P min )=2aL , (7) 

indicated by the vertical double-arrow in Fig.QJa). Here, 
fmin is the value of P(A\T, H) at the minimum "between 
the peaks" , while P max denotes the average peak value. 
The physical motivation for Eq. ([7]) is that, for cross- 
sectional areas "between the peaks", the bilayer reveals 
a coexistence between two slab domains where the to- 
tal interface length equals 2L [Fig. g^b)]. For DSPC, 
and assuming the lattice constant to be 1 nm, we ob- 
tain a ~ 1.1 pN, which is compatible with experimental 
values [11]. 

Next, we consider a DSPC bilayer with a fraction 
p = 0.03 of the lattice sites marked as "pinning sites". In 
Fig. Eta), we show distributions P(A\T, U) for T = 291 K 
obtained for 4 different random positions (samples) of 
pinning sites. Even though the temperature is the same 
as in Fig. 0Ja), a unique double-peaked distribution 
can no longer be identified. In contrast, P(A\T,H) is 
strongly sample dependent, and a multitude of rather 
exotic shapes is revealed. This behavior is characteris- 
tic of systems that belong to the universality class of the 
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FIG. 5. DSPC simulation results for T — 291 K in the 
presence quenched disorder, with a fraction of pinning sites 
p — 0.03. (a) Typical distributions P(A) for 4 different sam- 
ples of pinning sites, and system size L — 50. In contrast 
to Fig. [4[a) , a first-order transition can no longer be identi- 
fied, (b) Typical bilayer snapshot obtained at A = 54.4 A 2 
for a 200 x 300 lattice. A stable structure of multi-domains is 
observed. 



2D random- field Ising model (2D-RFIM) [2l[. Hence, 
by introducing the pinning sites, we have changed the 
universality class of the Pink model from ordinary 2D 
Ising toward 2D-RFIM (the pinning sites essentially cor- 
respond to a field of infinite strength acting at random 
locations) . 

There arc two features of the 2D-RFIM universality 
class that are remarkably consistent with experimental 
results for the main transition in supported membranes. 
First of all, 2D-RFIM universality implies the absence 
of macroscopic coexistence between gel and fluid do- 
mains [25l428l |. Indeed, inspection of simulation snap- 
shots [Fig. EJb)] reveals an equilibrium multi-domain 
structure, that is highly anisotropic, strongly resembling 
experimental AFM images [H(|. A second (related) fea- 
ture is that the 2D-RFIM has no true phase transition in 
the thermodynamic limit. In finite systems, there may 
be signs of a transition (or even several transitions; note 
that some of the distributions in Fig. [SJa) are triple- 
peaked) , but they will be "smeared" over a wide temper- 
ature range, and do not persist in the thermodynamic 
limit. Precisely this behavior has also been reported in 
experiments [4a , |50| . Simulation evidence that the pin- 
ning sites prevent a sharp transition in the thermody- 
namic limit follows from the (disorder-averaged) Binder 
cumulant Ui = [(A 2 )]/[(|A|) ], where [•] denotes an av- 
erage over different samples of pinning sites. As shown 
in Fig. [51 U\ —> 7r/2 as L increases, consistent with only 
a single phase. 



V. CONCLUSION 




293 T(K) 295 



FIG. 6. The (disorder-averaged) Binder cumulant U\ as a 
function of temperature T for DSPC with a fraction p — 0.03 
of pinning sites. In contrast to Fig. [2jb), an intersection of 
the curves for different L can no longer be identified. Instead, 
as the system size L increases, Ui —¥ n/2, indicative of the 
one phase region. For each system size, the disorder average 
comprised 200 samples of pinning sites. 



In this paper, the main phase transition in single- 
component phospholipid membranes was investigated us- 
ing the Pink model. Our simulations of the pure mem- 
brane (i.e. without quenched disorder) confirm the forma- 
tion of macroscopic gel and fluid domains below a critical 
temperature T c , and at the coexistence pressure IIcoex- 
We also demonstrated that, using the accepted values of 
the Pink model parameters, T c falls below experimentally 
measured transition temperatures. This explains why no 
phase transition was detected at the experimental transi- 
tion temperature in the simulations of Ref . [l8|. To resolve 
this issue, we propose that the strength of the van der 
Waals coupling in the Pink model be increased. By using 
the values proposed in this work, T c of the Pink model in 
the thermodynamic limit coincides with the experimen- 
tal main transition temperature. In addition, finite-size 
scaling was applied to confirm the universality class of 
the critical point in the Pink model, which was shown to 
be 2D Ising. This result holds irrespective of the number 
of conformational states q (as long as q > 2, of course). 
Hence, to capture the generic features of membrane phase 
behavior, a highly-detailed model is not always needed 
(which is consistent with the findings of Ref. |24T ) . 

We have also used the Pink model to describe the 
main transition in the presence of quenched disorder, 



8 



which may arise in case the membrane is deposited on a 
rough support. Assuming that this induces regions in the 
membrane where certain tail conformations become pre- 
ferred, the universality class changes toward that of the 
2D random-field Ising model. In the presence of quenched 
disorder, the Pink model reveals a stable multi-domain 
structure, and the absence of a sharp transition; these 
findings are indeed consistent with some of the experi- 
mental observations. 

Although the Pink model (both the pure version and 
the one containing quenched disorder) seems well suited 
to describe the main transition, we wish to end with a 
warning. By using the Pink model, one inevitably casts 
the main transition into the Ising universality class. This 
may not be entirely appropriate, as the main transition 
is essentially a melting transition leading to the forma- 



tion of nematic chain order. A liquid-crystal model may 
therefore be more suitable, such as the Maier-Saupe ap- 
proach followed in Rcf. [54j. If, indeed, the main transi- 
tion occurs close to a critical point [43|, l55| it may well be 
necessary to replace the discrete set of states of the Pink 
model by a continuous one In that case, one enters 
the regime of Hciscnbcrg-type models (which, provided 
certain conditions are met, do support 2D phase transi- 
tions p57l - [59j ) . The investigation of the main transition in 
terms of such a continuous model could be an interesting 
topic for future work. 
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